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Abstract 

We compare the Gram- Schmidt and covariant phase-space-basis-vector descriptions for three 
time-reversible harmonic oscillator problems, in two, three, and four phase-space dimensions re- 
spectively. The two-dimensional problem can be solved analytically. The three-dimensional and 
four-dimensional problems studied here are simultaneously chaotic, time-reversible, and dissipa- 
tive. Our treatment is intended to be pedagogical, for July 2011 publication in Communications 
in Nonlinear Science and Numerical Simulation, and for use in the second edition of our book on 
Time Reversibility, Computer Simulation, and Chaos. Comments are very welcome. 
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I. INTRODUCTION 

It is Lyapunov instability which makes statistical mechanics possible [1]. For stationary 
boundary conditions, either equilibrium or nonequilibrium, the exponential growth, oc e xt , 
of small phase-space perturbations {Sq, Sp} or {5q,5p,5(}, "sensitive dependence on initial 
conditions" , provides longtime averages independent of initial conditions. In addition to the 
coordinates and momenta { q,p }, time-reversible "thermostat variables" { C } can be used 
to impose nonequilibrium boundary conditions [2J, such as velocity or temperature gradients. 
A prototypical nonequilibrium problem simulates heat flow between two thermal reservoirs 
maintained at temperatures T# and Tq- The resulting heat flow through an internal New- 
tonian region, bounded by the two reservoirs, can then be studied and characterized 3, J]. 

To describe the Lyapunov instability for any of such systems, equilibrium or nonequi- 
librium, imagine the deformation of a small phase-space hypersphere <g>, comoving with, 
and centered on, a deterministic "reference trajectory" { q(t),p(t),((t) }• As the motion 
progresses the hypersphere will deform, at first becoming a rotating hyperellipsoid with 
the long-time-averaged exponential growth and decay rates of the principal axes defining 
the Lyapunov spectrum { A }. Much later the highly-distorted volume element, through 
the repeated nonlinear bending and folding caricatured by Smale's "Horseshoe" mapping, 
combined with an overall dissipative shrinking, comes to occupy a multifractal strange at- 
tractor. Though the apparent dimensionality of this steady-state attractor varies from point 
to point, it can be characterized by an overall averaged "information dimension" which is 
necessarily smaller, for stability of the phase-space flow (with ( (d/dt) In® ) < 0) than is 
the phase-space dimension itself [l, 5|. 

The time-averaged exponential growth and decay of phase volume have long been de- 
scribed by a set of orthonormal Lyapunov vectors { 5 }, one for each Lyapunov exponent 
A. Lyapunov spectra, sets of time- averaged local exponents, { A = ( A(i) ) }, for a variety 
of both small and large systems have been determined based on work pioneered by Stod- 

n n n 

dard and Ford|6], Shimada and Nagashima[7j, and Benettin's group[8i]. The summed-up 
Lyapunov spectrum can have thermodynamic significance, corresponding to the loss rate of 
Gibbs' entropy (the negative of the entropy gain of the thermal reservoir regions) in nonequi- 
librium systems interacting with Nose-Hoover thermostats, = S/k. The relative sizes 
of the phase-space components of the vector 5i associated with the largest instantaneous 
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Lyapunov exponent, Xi(t), allows instability sources ("hot spots", or better, "regions") to 
be located spatially. This localization of instability was long studied by Lorenz in his efforts 
to understand the predictability of weather. 

Recently ideas which had been expressed much earlier by Lorenz {9- 11] and Eckmann 
and Ruelle [l^| have been developed into several algorithms describing the phase-space de- 
formation with an alternative set of "covariant vectors" , vectors which "follow the motion" 
in a precisely time- reversible but somewhat arbitrary way 1 3M 1 91] . The literature describ- 



ing this development is becoming widespread while remaining, for the most part, overly 
mathematical (lots of linear matrix algebra) and accordingly hard to read. In order better 
to understand this work, we apply both the older Gram-Schmidt and the newer covariant 
time-reversibility ideas to three simple harmonic-oscillator problems. 

To begin, we describe the three example problems in Section II, along with the usual 
Gram-Schmidt method for finding local Lyapunov exponents and corresponding vectors. 
Some of the newer covariant approaches are outlined in Section III. Numerical results for 
the three problems, followed by our conclusions, make up the last two Sections, IV and V. 



II. LYAPUNOV SPECTRUM USING LAGRANGE MULTIPLIERS 
A. Model 1: Simple Harmonic Oscillator 

The models studied here are all harmonic oscillators. All of them incorporate variations 
on the textbook oscillator problem with coordinates, momentum p, and motion equations 
{ q = +p ; p = —q }• The first and simplest model[20| can be analyzed analytically. A (q,p) 
phase-space orbit of this oscillator is shown in Figure 1. Notice that the oscillator orbit 
shown there includes an arbitrary scale factor s, here chosen equal to 2. The corresponding 
Hamiltonian is H{q,p): 

2U = s +2 q 2 + s~ 2 p 2 ^ { q = +{p/A) ; p = -(4q) } ^q=-q. 

The "dynamical matrix" D for this oscillator describes the evolution of the set of infinitesimal 
offset vectors { 5 } needed to define the Lyapunov exponents, 

{ 5 = (g,p) sate llite - (<?, preference } ! { 5 = D ■ 5 } . 
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FIG. 1: Two views of a harmonic oscillator orbit with comoving ellipses centered on the orbit. 
The scaled equations of motion are q = +ps~ 2 ; p = —qs +2 . The changing aspect ratio of the 
ellipse shown at the left provides nonzero Gram-Schmidt Lyapunov exponents for the oscillator. 
The scaled plot at the right, of exactly the same data, shows that the exponents are a consequence 
of the scale factor s = 2 discussed in the text. 



D = 



In this case the matrix is 

' (1/4) 
-4 

Two local Lyapunov exponents, which reflect the short-term tendency of two orthonormal 
offset vectors to grow or shrink, can then be defined by the two vector relations: 

S 1 = D ■ Si - XnSi ; 

5 2 = D ■ 5 2 — X21S1 — X22S2 ■ 

We choose the "infinitesimal" length of the vectors, arbitrary in this linear problem, equal 
to unity for convenience. Snapshots of the two orthogonal Gram-Schmidt vectors are shown 
at the left in Figure 2. 

These local exponents have time-averaged values of zero, but sizable nonzero time- 
dependent and s-dependent fluctuations, 

( (A« - ( X tl » 2 ) = - s- 1 ) 2 ^ . 

The (scalar) Lagrange multipliers An and A 2 2 constrain the lengths of the offset vectors Si and 
S 2 while A 2 i constrains the angle between the vectors. The length-constraining multipliers 
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FIG. 2: Periodic phase-space orbit for the harmonic oscillator with scale factor s = 2. The Gram- 
Schmidt vector d{ is in the radial direction, identical to the first covariant vector 5f and necessarily 
perpendicular to 5 2 ■ The covariant vector 5% is parallel to the orbit. The Gram-Schmidt vectors 
are shown here at equally-spaced times and are, for this special model, identical in the two time 
directions, { 5$ } = { 5 b }. 

define the local Gram-Schmidt Lyapunov exponents. Their time averages give the global 
two Lyapunov exponents, Ai and A2 : 

A x = ( <Ji • D ■ 5 X ) = ( An ) ; A 2 = ( 5 2 ■ D ■ 5 2 ) = ( A 22 ) . 



This Lagrange-multiplier approach is the sma 
Gram-Schmidt approach to orthonormalization 21 



1 timestep limit of the finite-difference 
2 211. For a model with an iV-dimensional 



phase space iV Lagrange multipliers { Am } constrain the fixed lengths of the iV offset vectors 
while iV(iV — l)/2 additional multipliers { A^- } are required to keep the N(N — l)/2 pairs 
of vectors orthogonal. In the following two subsections we describe the analytic formulation 
of the Lagrange-multiplier problem for examples with three- and four- dimensional phase 
spaces. For large N the Gram-Schmidt approach is considerably faster and simpler than the 
Lagrange-multiplier one. 



B. Model 2: Harmonic Nose-Hoover Oscillator with a Temperature Gradient e 

The second model, unlike the first, can exhibit long-term chaotic motion, with three- 
dimensional phase-space perturbations { 8q, Sp, 5( } growing or shrinking exponentially in 
time, oc e xt . Here ( is a friction coefficient and controls the instantaneous changing kinetic 
temperature p 2 , so as to match a specified target temperature T(q) [231 124|. Here we allow 



the temperature to depend upon the oscillator coordinate [25l428(|. 

1 - e < T(q) = 1 + etanh(g) < 1 + e . 

The temperature gradient makes overall dissipation possible, characterized by a shrinking 
phase-space volume, (g) — > 0, and resulting in a strange attractor, with Dj < 3, or even a 
one-dimensional limit cycle, Dj = 1. 

The e-dependent temperature gradient opens the irresistable possibility for heat to be 
absorbed at higher temperatures than those where it is expelled. The Nose-Hoover equations 
of motion for the nonequilibrium oscillator are : 

{q = p; p = -q-(p; ( = p 2 - T(q) } . 

The friction coefficient ( allows the long-time-averaged kinetic temperature, proportional to 
( p 2 ), to conform to a nonequilibrium steady state, characterized by the imposed tempera- 
ture profile T(q) : 

( ( ) constant — > ( p 2 ) = ( T(q) ) . 

Because the motion occurs in a three-dimensional phase space the dynamical matrix D, 
which governs the motion of phase-space offset vectors is 3 x 3 : 

(dq/dq) (dq/dp) (dq/d()] \ 1 
D = (dp/dq) (dp /dp) {dp/dQ = —l—^—p 
(dC/dq) (dC/dp) (d(/d()\ [ -T 2p 
where T' is the derivative of temperature with respect to q : 

-(OC/dq) = V = ecosh- 2 (g) . 

For simplicity, we choose these vectors to have unit length. In addition, six Lagrange mul- 
tipliers are required to maintain the orthonormality of the three offset vectors { S±, 62,63 } 

A n 
A21 A22 
A31 A32 A33 

{ A x = ( \ u (t) ) ; A 2 = ( A 32 (f) ) ; A 3 = ( A 33 (t) ) } . 

This model can exhibit chaos or regular behavior, depending on the initial conditions as well 
as the maximum temperature gradient e. Chaotic solutions for this oscillator are necessarily 
numerical, rather than analytical, and are illustrated in Section III. 



C. Model 3: Doubly Thermostated Oscillator with a Temperature Gradient e 

The Nose- Hoover oscillator, with a single friction coefficient ( (Model 2), is not ergodic. 
For relatively small temperature gradients it exhibits an infinity of regular solutions bathed in 
a chaotic sea. For a glimpse of the details see Reference 24. This geometric three-dimensional 
complexity can be reduced at the price of introducing a second thermostat variable. The 
resulting four- dimensional model, as well as a similar extension treated in Reference 18, 
has two friction coefficients (CO rather than just one. With a constant temperature T the 
resulting ergodic canonical-ensemble phase-space probability density is 

f(q,P,(,0 = (l/^T)e-^ T e-^ T e-^e^ 2 / 2 . 

The nonequilibrium multifractal extension of this model results if (as in Model 2) the temper- 
ature depends upon the oscillator coordinate, T(q) = 1 + etanh(g). In this nonequilibrium 
case there are four equations of motion: 

{ q = P ; P = ~q - (P ~ & 3 ; C = P 2 - T (q) ; £ = P* - 3p 2 T(g) } . 

The corresponding four- variable dynamical matrix D, which controls the linearized 
( "tangent-space" ) motion, 5 — D • 5, is 

10 

D= -i [-C - ^P 2 ] ~P -p 3 
-V 2p 

-3p 2 T' [4p 3 -6pT] 

A lower-triangular array of constraining Lagrange Multipliers, can then be defined. Just 
as before, the diagonal multipliers maintain the lengths of the vectors constant and the 
off-diagonal multipliers maintain the orthogonality of the vectors. For this model, with four 
offset vectors, we have ten Lagrange multipliers in all: 

5 1 = D ■ 5i - \ n 5i ; 
8 2 = D ■ 62 — X21S1 — A 22 <5 2 ; 

63 = D ■ 63 — \si5i — A32$2 — ^33^3 ; 
$4 = D ■ $4 — A4i$i — A42$2 ~ A43$3 — A44$4 . 
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The Gram- Schmidt Lyapunov exponents are the long-time-averaged diagonal Lagrange 
Multipliers I22I, |23j : 



Aj = (A«(t)) = (St ■ D ■ Si) /{Si ■ 6,) = (Si-D- S t ) . 

Here, as is usual, we choose the (arbitrary) length of the tangent-space offset vectors equal 
to unity: { \S\ = f }. The instantaneous off-diagonal elements, 

{ \ ij (t) = 6 i -D-6 j + S r D-6 i } , 

describe the tendency of the offset vectors to rotate relative to one another. This four- 
dimensional model is already sufficiently complex to illustrate the distinctions between the 
Gram-Schmidt and covariant descriptions of phase-space instabilities. We include numerical 
results for this model too, in Section IV. 



III. COVARIANT LYAPUNOV VECTORS AND THEIR EXPONENTS 



Several loosely-related schemes for evaluating covariant (as opposed to Gram-Schmidt) 
Lyapunov vectors have been described and explored. Their proponents are mostly interested 
in the mathematics of weather modeling predictability or in better understanding the sta- 
tistical mechanical sensitivity to phase-space perturbations. Lorenz' famous 1972 talk title: 
"Predictability: Does the Flap of a Butterfly's Wings in Brazil Set off a Tornado in Texas?" 



indicates the broad scope of this work [111]. 

Wolfe's 2006 dissertation [ljl and the Kuptsov-Parlitz review|l9j are particularly useful 
guides to this work. Some of the underlying ideas date back to Lorenz' early studies, in 
1965 [9j. The main idea is not so different to the old Gram-Schmidt approach and in fact 
requires information based on that Gram-Schmidt (or Lagrange multiplier) approach not 
only forward (as is usual), but also backward (which can be artificial, violating the Second 
Law of Thermodynamics) in time. 

The main idea is to seek a representative basis set of comoving and corotating infinitesimal 
phase-space vectors { S c }, (c for covariant) guided by the linearized (treating { q,p,(i£ } 
as constants in the D matrix) motion equations, 

{ 5 C = D ■ 5 C } . 
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The covariant basis vectors follow the linearized flow equations, without Lagrange multipli- 
ers, and so are generally not orthogonal. Provided that the "unstable" manifold, made up 
of the phase-space directions corresponding to longtime expansion, can be usefully distin- 
guished from the "stable" manifold, corresponding to directions associated with longtime 
contraction, the covariant basis vectors are locally parallel to these two manifolds. 

A "covariant" set of vectors would seem not to require Gram-Schmidt constraints, or 
Lagrange multipliers, because orthogonality is not required. Even so, all the existing com- 
putational algorithms which have been developed to find covariant vectors begin by finding 
the orthonormal Gram-Schmidt vectors. Two sets of Gram-Schmidt vectors are the next 
requirement, the usual forward-in-time vectors { 8* } and the not-so-usual backward- in- 
time set { 5 b }. This second set of vectors is obtained by post-processing the time-reversed 
phase-space trajectory. 

In time- reversible energy-conserving Hamiltonian mechanics the reversed trajectory can 
be as "natural" as the forward trajectory. In dissipative systems, with phase-space shrinkage, 
the stored and reversed trajectory is typically "unnatural", and would violate the Second 
Law of Thermodynamics. Unlike the Gram-Schmidt forward and backward vectors [ we 
will denote them by { 5* } and { 5 b } ], the covariant vectors { 5 C } are defined so as to be 
identical, at least for Hamiltonian mechanics, in the two time directions, that is, "covariant" . 
The Gram-Schmidt behavior, with the forward and backward vectors generally different, is 
a symmetry breaking whose source is not at all apparent in the underlying time-symmetric 
differential equations of motion. 

The differential equations for the time development of the covariant vectors account 
for the simultaneous stretching and rotation of an infinitesimal phase-space hypersphere. 
Evidently the maximum growth (and maximum shrinkage) rates correspond to the usual 
largest Lyapunov exponents forward and backward in time, 8{ = 81 and 5\ = 5 C N in an 
N- dimensional phase space. Another covariant offset vector parallels the local direction of 
the phase-space velocity, (q,p, (). The remaining N — 3 covariant vectors require more work. 
The difference between the Gram-Schmidt and covariant phase-space vectors is illustrated 
for the simple harmonic oscillator in Figure 2. 

The relatively readable Wolfe-Samelson approach [if]] begins with the Gram-Schmidt sets 
{ 5^,5 b }, forward and backward in time. Then covariant unit vectors are associated with 
the long-time-averaged individual Lyapunov exponents, beginning with the second (or, using 
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time symmetry, beginning with the next-to-last and working backwards). The covariant 
vectors are then expressed as linear combinations of (some of) the forward and backward 
Gram-Schmidt vectors. 

The problem is overdetermined, in that 2N Gram-Schmidt vectors are used as bases for 
the N covariant vectors. The second covariant vector can be written as a linear combination 
of the first two Gram-Schmidt vectors: 

$2 = y(4 + yl4 . 

The constants { yf } are then determined by solving a relatively simple eigenvalue problem. 
Likewise, the next-to-last covariant vector can be written in terms of the first two Gram- 
Schmidt vectors from the time-reversed trajectory: 

*at-i = y b A + y^2 ■ 

In general the nth covariant vector can be expressed as a sum of n Gram-Schmidt vectors 
with the constants y determined by solving a set of linear equations. All of these vectors 
are unit vectors, with length 1. The many method variations (using the first few, the last 
few, or some of both sets of Gram-Schmidt vectors) are discussed in Wolfe's thesis, which is 
curre ntly available online. The Appendix of Romero-Bastida, Pazo, Lopez, and Rodriguez' 
work [17], as well as the Kuptsov-Parlitz review 19[ are also useful guides. 



The main difficulty in putting all of this work into perspective is a result of the frac- 
tal/singular nature of the phase-space vectors. This structure can be traced to bifurcations 
in the past history and/or in the future evolution of a particular phase point. This sensitivity 
to initial conditions means that slightly different differential-equation algorithms can lead 
to qualitatively different trajectories, making it hard to tell whether or not two computer 
programs are consistent with one another. The best approach to code validation is to repro- 
duce properties of the covariant vectors which are insensitive to the integration algorithm. 
We will consider some of these properties for our three example oscillator models. 
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IV. RESULTS FOR THE THREE HARMONIC OSCILLATOR PROBLEMS 



A. Ordinary One-Dimensional Harmonic Oscillator with Scale Factor s = 2 

The first of our three illustrative oscillator problems is an equilibrium problem in Hamil- 
tonian mechanics, a harmonic oscillator with unit frequency and with a scale factor s set 
equal to 2 : 

q=-q^m = s +2 q 2 + s~ 2 p 2 ^ { q = +(p/4) ; p = -(4g) } — ► q = -q . 

The local Lyapunov exponents for such an oscillator depend upon s and vanish in the usual 
case where s = 1 . 

A (q,p) phase-space plot of a typical orbit for s = 2 was shown in Figure 1. The oscillator 
orbit we choose to analyze is the ellipse shown there: 

{ q = cos(t) ; p = -4sin(t) } ; lQq 2 + p 2 = 16 = 2H . 

The local Lyapunov exponents describe the instantaneous growth rates of infinitesimal "off- 
set vectors" comoving with the flow. These are of three kinds, orthonormal vectors moving 
forward in time { 8* }, orthonormal vectors moving backward in time { 8 b }, and special co- 
variant vectors { 5 C } whose forward and backward orientations in phase space are identical. 
The initial direction of each of these vectors is arbitrary because a reversed (antiparallel) 
vector (+8 i — > —8) satisfies exactly the same equations and definitions. 

For the simple scaled harmonic oscillator it is possible to solve analytically for the Gram- 
Schmidt vectors forward and backward in time as well as the covariant vectors and all the 
corresponding local exponents, 

{ 8f(t),8\t),8 c (t) } — > { \f(t),\\t),\ c (t) } , 

once the initial conditions are specified. See Figure 2 for an illustration of the following 
choice of initial values, at t — : 

(q,p) = (1,0) ; S{ = 8\ = 8{ = (1,0) ; 8 f 2 = 8 b 2 = 8 C 2 = (0, 1) . 

In the course of the motion, with period 2n, the first forward Gram-Schmidt vector S(, which 
is generally identical to the first covariant vector, turns out to be always radial. This vector 
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is necessarily perpendicular to 82(f). The vectors 8[ and 8% generally differ. The second 
covariant vector is a unit vector which remains always parallel or antiparallel (and hence 
"covariant") to the orbit, 

8 c 2 = ±(q,p)/\(q,p)\ . 

Generally, all the covariant vectors satisfy the linearized small-5 motion equations. Because 
the oscillator motion equations are linear, these oscillator results are correct for finite, not 
just infinitesimal, vectors. 

The Gram-Schmidt Lyapunov exponents can be calculated as Lagrange multipliers which 
maintain the orthonormal relations between { 8\(t), 82(f) } in both the forward and the 
backward directions. For this simple problem the vectors are unchanged by time reversal, 
S{ = 8\ ; 8{ = 8\. The first vector is not constrained in orientation, but is restricted to 
constant length by the Lagrange multiplier An: 

{ Sqi = (l/4)5pi - \n8qi ; 8 Pl = -ASqi - XuSpi } — > An = -(15/A)Sq 1 Sp 1 . 

Using An m either the 8q\ or 8p\ equation gives a differential equation for the unit vector 
81 = (8qi,8pi), which simplifies with the definition of a new variable, the angle 9: 

5qi = (V4)5pi + (lh/A)8 Pl 8ql 5 S Pl = -4<J 9l + (15/4)^^ 

[ 8q\ = cos(9) ; 8pi = — sin(6 l ) ] — > 

= (1/4) + (15/4) cos 2 (^) = 4 - (15/4) sin 2 (^) — > 9 = arctan(4tan(t)) . 

Although the probability density for 8q [ or 8p ] diverges, whenever q = [ or whenever 
p — ], the probability density for 9, which describes the nonuniform phase-space rotary 
motion, is well-behaved, as shown in Figure 3. 

The instantaneous Lyapunov vectors, both Gram-Schmidt and covariant, from the orbit 
shown in Figure 2, are shown as functions of time in Figure 4. The exponents are: 

A^ (t) = \ b (t) = ±(-15/4) sin(0) cos(#) ; 9 = arctan(4tan(t)) , 

with the second covariant exponent given by the strain rate parallel to the orbit: 

dhiv\\/dt = -15sin(t)cos(t)/[sin 2 (t) + 16cos 2 (t)] . 

Figure 4 shows that this covariant exponent is 90° out of phase with the exponent associated 
with 8((t). 
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Figure 3 
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FIG. 3: Probability density for the angle 9 = arctan(4tan(t)) for the scaled harmonic oscillator. 

Figure 4 
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FIG. 4: Lyapunov exponents for the harmonic oscillator with scale factor 2. The Gram-Schmidt 
exponents are the two heavy lines, and sum to zero. The second covariant exponent, A^t) = 
(d In /dt) along the orbit, is the lighter dashed line. 

Because the oscillator is not chaotic, these detailed results depend upon the choice of 
initial conditions. The general features illustrated by this problem include (1) the orthogo- 
nality of the Gram-Schmidt vectors, (2) the possibility of obtaining a "reversed" trajectory 
by using a stored forward trajectory, changing only the sign of the timestep +dt — > —dt, 
(3) the identity of the first Gram-Schmidt vector with the first covariant vector, and (4) 
the identity of the trajectory direction with the covariant vector corresponding to the time- 
averaged exponent ( An ) = 0. In the next problem we study the directions of the forward and 
backward vectors and show that they typically differ, though the time- averaged exponents 
(for a sufficiently long calculation) agree, apart from their signs. 
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Figure 5 




FIG. 5: Phase-space trajectories for a Nose-Hoover heat conducting oscillator with the coordinate- 
dependent temperature, T = 1 + etanh(g). The strange attractor at the left corresponds to a 
maximum temperature gradient e = 0.25 while the limit cycle at the right corresponds to e = 0.50. 
Here the coordinate q increases from left to right, the momentum p from front to back, and the 
friction coefficient ( from bottom to top. 

B. Three-Dimensional (q,p,() Nose-Hoover Oscillator 

Here we consider the second model oscillator. Its nonequilibrium temperature profile, 

T(q) = 1 + etanh(g) , 

provides both chaotic and regular solutions, depending on the strength of the temperature 
gradient. Trajectory segments for two values of e (the maximum value of the temperature 
gradient) are shown in Figure 5. For the smaller value e = 0.25, the second Lyapunov 
exponent vanishes, A2 = ( A^t) ) = 0. Its probability density, shown in Figure 6, is quite 
different to that of the second covariant exponent, A^t), which remains parallel to the 
trajectory direction at all times. For the larger limit-cycle value, e = 0.50, the longtime 
orbit is a limit cycle, with period 8.650. On this limit cycle the largest Lyapunov exponent 
is identical to the largest covariant exponent, and necessarily vanishes. The time dependence 
of this exponent, \{i(t) = A'[ 1 (t), with 

(\{ 1 (t)) = (\l 1 (t)) = \{ = \l = 0, 

is shown in Figure 7. 

The three-dimensional system of motion equations, 

{ q = p ; p = -q-Cp; C = p 2 - T (q) } , 
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FIG. 6: Logarithms of probability densities, over about six decades, for the orbital strain rate, 
^22 (*) = (dlnv/dt) and the corresponding instantaneous Gram-Schmidt Lyapunov exponent, 
for the chaotic oscillator with e = 0.25. 

Figure 7 
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FIG. 7: Instantaneous value of the orbital strain rate for the limit cycle with e = 0.50. The strain 
rate is identical to the corresponding instantaneous Gram-Schmidt Lyapunov exponent, X 1± . 

can be "reversed" in either of two ways: (1) replace +dt by — dt in the fourth-order Runge- 
Kutta integrator or (2) leaving dt > unchanged, change the signs of the momentum p and 
the friction coefficient ( so that the underlying physical system traces its coordinate history 
backward in time, q(+t) — > q(—t). In either case longtime stability with decreasing time 
requires that the dissipative forward trajectory be stored and reused. If the trajectory is not 
stored then the reversed motion abruptly leaves the unstable reversed orbit. Figure 8 shows 
the jump from the illegal reversed motion, violating the Second Law of Thermodynamics, 
to a stable obedient motion after sixteen circuits of the illegal cycle. 

For simplicity of notation and description we have adopted the first reversibility definition 
throughout the present work. It is necessary to recognize that the "motion" (forward in time) 
we analyze is dissipative, though time-reversible, while the reversed motion (backward in 
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Figure 8 




FIG. 8: Time history of a Nose-Hoover oscillator in a temperature gradient, e = 0.50. The 
variation of { p, ( } in time just after the reversal at time zero is shown. After 16 "wrongway" 
limit cycles, during which entropy is absorbed by the oscillator, the thermodynamically unstable 
trajectory appears to jump to a cycle which obeys the Second Law of Thermodynamics. The jump 
time of 150 closely corresponds to the estimate appropriate to 48-bit precision, a cycle period of 9, 
and the known Lyapunov exponent of 0.222: e xt = e 222 ' = 2 48 — > t = 150 ~ 16 x 9. 

time and contrary to the Second Law of Thermodynamics) is actually so unstable ( because 
( <8> ) > ) as to be unobservable unless the trajectory has been stored in advance. 

The instability of the reversed trajectory, leading to symmetry breaking, can be under- 
stood by considering the flow of probability density in phase space. The comoving proba- 
bility density f(q,p, (, t) obeys the analog of Liouville's continuity equation. At any instant 
of time the changing phase-space probability density responds to the friction coefficient ( : 

[///] = (dq/dq) + (dp/dp) + (OC/dC) = + C + = -®/<g> = -X f n (t) - \{ 2 {t) - \Ut) . 

Phase- volume expands/contracts when ( is negative/positive. The time- averaged value of 
the friction coefficient ( is necessarily positive, and reflects the heat transfer (from larger to 
smaller values of the oscillator coordinate q) consistent with the Second Law of Thermody- 
namics. 

Because this Nose-Hoover problem is chaotic, no analytic solution is available, making 
it necessary to compute the vectors and exponents numerically. Lyapunov exponents have 
long been determined numerically. Spotswood Stoddard and Joseph Ford's pioneering im- 
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plementation was generalized by Shimada, Nagashima, and by Bennetin's group. These 
latter authors used Gram-Schmidt rescaling of phase-space offset vectors. This continuous 
limit of this numerical approach was formalized by introducing Lagrange multipliers { \j } 



to constrain the offset vectors' orthonormality 
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This problem raises a computational question: How best to distinguish the usual Gram- 
Schmidt vectors { 8* } from the covariant ones { 5 C }? Some quantitative criterion is 
required, at the least for reproducibility if not for clarity. In the case of the thermostated 
oscillator we have computed the mean-squared projections into (q,p, C) space for the three 
Gram-Schmidt basis vectors, both forward and backward in time. The results are essentially 
the same for either time direction: 

(0.22, 0.29, 0.49), (0.27, 0.38, 0.35), (0.51, 0.32, 0.17) for e = 0.25 ; 

(0.22, 0.25, 0.53), (0.27, 0.41, 0.32), (0.51, 0.34, 0.15) for e = 0.50 . 

Evidently, at least for this problem, there is little geometric dependence of the Gram-Schmidt 
vectors on the presence or absence of chaos. On the other hand, the second covariant 
vector (parallel to the flow) has mean-squared components (0.21, 0.28, 0.51) for e = 0.25 
and (0.22, 0.25, 0.53) for e = 0.50, and so is significantly different to the Gram-Schmidt 
#2 or 8\. A shortcoming of work on the covariant vectors so far has been the lack of 
quantitative information concerning them. This is an undesirable situation, as it makes 
checking computations problematic. 

We turn next to a four- dimensional problem. Although the first and last Gram-Schmidt 
vectors and the trajectory direction give three of the four covariant vectors, the fourth 
requires new ideas, and illustrates the general case of determining a complete set of covariant 
basis vectors. 



C. Four-Dimensional (q,p, CO Doubly-Thermostated Oscillator 

Problems with four or more phase-space dimensions require the mathematics, mostly lin- 
ear algebra, of covariant vectors for their solution. Here we choose a doubly-thermostated 
harmonic oscillator (two thermostat variables, ( and £), again with a nonequilibrium tem- 
perature profile, T = 1 + tanh(g). The equations of motion (which give the canonical 
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phase-space distribution characteristic of the temperature T when e vanishes) are: 

{ q = v ; v = -q - Cp - & 3 ; C = p 2 -T ; £ = p 4 - 3p 2 T } . 



Again T is the kinetic temperature, the time-averaged value of p 2 . ( and £ are time- 
reversible thermostat variables which control the second and fourth moments of the velocity 
distribution. At equilibrium the solution of these motion equations is the complete canonical 
phase-space distribution. For T = 1 this stationary distribution has the form: 

f(q, p, C, oc exp[-(g 2 + p 2 + ( 2 + £ 2 )/2] ■ 

This oscillator system is ergodic. When the temperature T is made to depend upon the 
coordinate q, a nonequilibrium strange attractor results. We can give some detailed results 



for the well-studied 



25H27I] special case < T(q) = l + tanh(g) < 2. For this nonequilibrium 



problem the Lyapunov spectrum forward in time is already known, 

{ \f } = { +0.073, 0.000, -0.091, -0.411 } — > { \ b } = { +0.411, +0.091, 0.000, -0.073 } , 

as is also the strange attractor's phase-space information dimension, Dj = 2.56, reduced by 
1.44 from the equilibrium case (e = 0), where the spectrum is symmetric: 

{A} = { +0.066, 0.000, 0.000, -0.066 } . 



This information dimension result is in definite violation of the Kaplan- Yorke conjecture 271 ]. 

Di = 2.56 = D KY = 2 + (0.073/0.091) = 2.80 . 

The Gram-Schmidt vectors { 8*, 5 b } for this problem are easily obtained using the La- 
grange Multiplier approach, and give three of the four covariant vectors: 



51 = 51; 8 c 2 = v/\v\; 5 C A = 5\ 



b 



where v is the phase-space trajectory velocity, (q,p, C, £)• 

Wolfe and Samelson show that the only missing covariant vector, 5|, can be expressed as 
a linear combination of the first three forward (or the last two backward) 5 vectors. Either 
choice requires solving an eigenvalue problem numerically and both choices lead to exactly 
the same covariant vector 8%, even if the Gram-Schmidt vectors are not perfectly converged. 
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Figure 9 




FIG. 9: Probability densities for A^i) and A^i)- The logarithmic ordinate scale covers two 
decades with data from the central two billion out of four billion fourth-order Runge-Kutta 
timesteps of 0.001 each. 

A comparison of the various local Lyapunov exponents is detailed in Figures 9-12. We 
have ordered the thickness of the plotted lines throughout, in the order 1>2>3>4so 
that the thickest line represents the distribution for \{i(t) (mean value +0.073) and, in the 
Backward plot, \\i{t) (mean value +0.411). Visually, the data for two billion points, shown 
here, are indistinguishable from data for two hundred million points, shown in Version 1 of 
this work. 

With a single exception, the vanishing Gram-Schmidt exponent derived from S^t), the 
data suggest a fractal nature for the probability densities of the various exponents. Figure 
9 compares the probability distributions for A^(^) and the covariant exponent A^^) corre- 
sponding to motion along the trajectory. The Gram-Schmidt histogram is much smoother 
than that of its covariant cousin. At the moment we have no explanation for this interesting 
qualitative difference. 

Otherwise the local exponent distributions are qualitatively similar — large fluctuations 
compared to the actual exponent values. The extra work associated with the local covariant 
spectrum cannot be justified based on these data. It should be noted that the Gram-Schmidt 
exponents forward and backward in time are quite different, reflecting the difference between 
the future and the past trajectories. 

To make it possible for the reader to coordinate our results with those of Wolfe and 



Samelson 



we write a set of equations which can be solved in order to find the third 
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covariant vector 5% as an expansion in the forward-in-time Gram-Schmidt vectors: 

nti-wi-so nrt-s b k M-sb nti-wiH)' 

U4-s b k M.5[) u4-s b k M-4) u4-5 b k )(5 b k -4) ■ 

nu-m b k-H) nu-miti) nu-%M-U)_ 

Here the sums include the two values of k — 1, 2. The numbering system, though arbitrary, 
is crucial. Forward in time the correspondence of the vectors with the long-time-averaged 
Lyapunov exponents is 

5( ->■ +0.073 ; 8{ ->■ 0.000 ; 5{ ->■ -0.091 ; 5{ ->■ -0.411 . 

Backward in time the ordering is the same with the signs changed: 

5\ ->■ -0.073 ; ^ ->■ 0.000 ; <?| ->■ +0.091 ; $ ->■ +0.411 . 

Of course the £/&mZ covariant vector going forward could equally well be viewed as the 
second going backward using an expansion in terms of the backward-in-time Gram-Schmidt 
vectors. This alternative approach results in a smaller matrix without the sum over k: 

(51 ■ 5()(5( ■ 51) (51 ■ 5()(5( ■ 5 b 2 ) ' 
(5\ ■ 5{)(5{ ■ 51) (5\ ■ 5{)(5{ ■ 5 b 2 ) _ 

Here the ordering of the Gram-Schmidt vectors follows this same new ordering of the back- 
ward vectors: 

51 ->■ +0.411 ; 5 b 2 ->■ +0.091 ; 5 b 3 ->■ 0.000 ; 8\ ->■ -0.073 . 

The vectors forward in time follow the same ordering with the signs changed. 

Once the ordering of the vectors is correctly negotiated one can find the eight Gram- 
Schmidt vectors (four in each time direction) and the four covariant vectors (the same in 
either time direction). For comparison we show the probability densities for the vectors in 
Figures 10, 11, and 12. 

V. PERSPECTIVES 

The covariant phase-space vectors, illustrated here for oscillators, are said to have two 
advantages: [1] they display the time-symmetry of the underlying equations of motion and 
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M f ■ y = ; M f = 



M b ■ y = ; M b = 



Figure 10 







^^^^^ward Gram-Schmidt Lyapunov Exponijnts^^! 







-2 < X < +2 

FIG. 10: Probability density for the four forward-in-time Lyapunov exponents, { A^(t) }. Notice 
that the distribution for the largest of the exponents (widest line) is the same as that for the 
covariant Lyapunov exponent Af 1 (i). The logarithmic ordinate scale covers two decades with data 
from the central two billion out of four billion fourth-order Runge-Kutta timesteps of 0.001 each. 

Figure 11 




-2 < X < +2 



FIG. 11: Probability density for the four backward-in-time Lyapunov exponents, { A^(t) }. Notice 
that the probability distribution for the most positive of the backward exponents, \\ = +0.411 
(widest line) matches that for the covariant Lyapunov exponent A| 4 (i) shown in Figure 12. The 
logarithmic ordinate scale covers two decades with data from the central two billion out of four 
billion fourth-order Runge-Kutta timesteps of 0.001 each. 

[2] they provide results which are "norm-independent". The Gram-Schmidt vectors differ, 
reflecting the past and steadfastly ignorant of the future. Consider a purely-Hamiltonian 
situation, the sudden inelastic collision of two rapidly-moving blocks, converting ordered 
kinetic energy to heat. The Gram-Schmidt vectors turn out to be more localized in the 
forward direction of time than in the (completely unphysical) reversed direction, where 
the particle trajectories have been stored [29]. The phase-space directions corresponding 
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Figure 12 















Covariant Lyapunov Exponents"— 


-2 < 


X 


< +2 



FIG. 12: Probability density for the four covariant Lyapunov exponents, { A^(t) }. The largest 
and smallest of these exponents match the largest forward and the largest backward of the Gram- 
Schmidt exponents. The covariant vectors are identical forward or backward in time, with their 
exponents simply changing sign. The logarithmic ordinate covers two decades with data from the 
central two billion out of four billion fourth-order Runge-Kutta timesteps of 0.001 each. 

to maximum growth and decay (easy to calculate by singular value decomposition of the 
symmetrized D matrix and relatively smooth in phase space, rather than multifractal) show 
the same time symmetry as do the covariant vectors. This time symmetry seems unhelpful 
for describing the irreversibility of this purely-Hamiltonian shockwave process. We like the 
time asymmetry of the Gram-Schmidt vectors, reflecting as they do the past-based nature 
of the Second Law of Thermodynamics. 

Despite the norm-independence of the Lyapunov vectors it is apparent that "local" (in 
phase-space or in time) Lyapunov exponents depend on the chosen coordinate system. Even 
the simple harmonic oscillator considered here shows coordinate-dependent local exponents. 

The time symmetry-breaking associated with time-reversible dissipative systems (like the 
three-dimensional and four- dimensional oscillator models studied here) prevents backward 
analyses unless trajectories are stored. Lyapunov instability helps attract a repellor state 
(violating the Second Law) to the Law-abiding strange attractor. Like the overall trajec- 
tory, the Gram-Schmidt vectors are themselves unstable to symmetry breaking. A reversed 
trajectory calculation will not also reverse the Gram-Schmidt vectors unless those vectors 
too are stored from a forward trajectory. The covariant vectors are the same going forward 
and backward in time. The symmetry breaking which naturally occurs can only be avoided 
by repeating a trajectory, making it possible to find the covariant solution. 



22 



Despite the double-precision accuracy of Runge-Kutta solutions it is difficult to validate 
computer programs for dissipative chaotic systems precisely due to their Lyapunov instabil- 
ity. For this reason the squared projections of the offset vectors and the exponent histograms 
are good choices for validation. The details of specific program results necessarily vary and 
reflect the fractal nature of the phase-space singularities describing the inevitable past and 
future bifurcations. The fractal nature of these singular points frustrates any attempt to 
gain accuracy through mesh refinement. Both the covariant and the Gram-Schmidt expo- 
nents share this fractal nature, but differently. The Gram-Schmidt vectors are analogous 
to passengers' reactions to a curvy highway while the covariant view is that of a stationary 
pedestrian observer. 
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